#include <iostream>
#include <cmath>
#include "CmdLine/CmdLine.hh"
#include "scalar.h"

#include "TFile.h"

using namespace std;

int main(int argc, char** argv)
{
  CmdLine cmdlin(argc, argv);

  //200 GeV scalar
  Scalar mymodel(500, 0, 1);

  //find minima
  for(double T=50; T<500; T+=50)
    {
      cout<<"(T, vev): "<<T<<", "<<mymodel.find_v(T)<<endl;
    }


  double vev=0;

  double Tc=mymodel.Tc(vev);
  cout<<"Tc: "<<Tc<<", "<<vev<<endl;

  //now scan and find all the Tc
  vector<double> x1;
  vector<double> x2;
  vector<double> x3;
  vector<double> y1;
  vector<double> y2;
  vector<double> y3;

  for(double s_mass=80; s_mass<400; s_mass+=10)
    {
      //declare a new model
      Scalar temp_model(s_mass, 0, 1);
      
      //get the critical temperature
      double vevp=0;
      double Tp=temp_model.Tc(vevp);

      double order=0;
      if(Tp>0) order=vevp/Tp;

      if(order >0)
	{
	  x1.push_back(s_mass);
	  y1.push_back(order);
	}
    }

  for(double s_mass=80; s_mass<400; s_mass+=10)
    {
      //declare a new model
      Scalar temp_model(s_mass, 0, 2);
      
      //get the critical temperature
      double vevp=0;
      double Tp=temp_model.Tc(vevp);

      double order=0;
      if(Tp>0) order=vevp/Tp;

      if(order >0)
	{
	  x2.push_back(s_mass);
	  y2.push_back(order);
	}
    }

  for(double s_mass=80; s_mass<400; s_mass+=10)
    {
      //declare a new model
      Scalar temp_model(s_mass, 0, 3);
      
      //get the critical temperature
      double vevp=0;
      double Tp=temp_model.Tc(vevp);

      double order=0;
      if(Tp>0) order=vevp/Tp;

      if(order >0)
	{
	  x3.push_back(s_mass);
	  y3.push_back(order);
	}
    }


  TFile fout("out.root", "RECREATE");

  TGraph* phase=new TGraph(x1.size(), &x1[0], &y1[0]);
  phase->SetName("phase1");
  phase->Write();

  TGraph* phase2=new TGraph(x2.size(), &x2[0], &y2[0]);
  phase2->SetName("phase2");
  phase2->Write();

  TGraph* phase3=new TGraph(x3.size(), &x3[0], &y3[0]);
  phase3->SetName("phase3");
  phase3->Write();

  TGraph* gc=mymodel.GetGraph(Tc);
  TGraph* g0=mymodel.GetGraph();
  TGraph* g20=mymodel.GetGraph(20);
  TGraph* g40=mymodel.GetGraph(40);
  TGraph* g60=mymodel.GetGraph(60);
  TGraph* g80=mymodel.GetGraph(80);
  TGraph* g100=mymodel.GetGraph(100);
  TGraph* g120=mymodel.GetGraph(120);
  TGraph* g140=mymodel.GetGraph(140);
  TGraph* g160=mymodel.GetGraph(160);
  TGraph* g180=mymodel.GetGraph(180);
  TGraph* g200=mymodel.GetGraph(200);
  TGraph* g250=mymodel.GetGraph(250);
  /*
  TGraph* g300=mymodel.GetGraph(300);
  TGraph* g400=mymodel.GetGraph(400);
  TGraph* g500=mymodel.GetGraph(500);
  TGraph* g600=mymodel.GetGraph(600);
  TGraph* g700=mymodel.GetGraph(700);
  TGraph* g800=mymodel.GetGraph(800);
  TGraph* g900=mymodel.GetGraph(900);
  TGraph* g1000=mymodel.GetGraph(1000);
  */

  gc->SetName("gc");
  gc->Write();
  g0->Write();
  g20->Write();
  g40->Write();
  g60->Write();
  g80->Write();
  g100->Write();
  g120->Write();
  g140->Write();
  g160->Write();
  g180->Write();
  g200->Write();
  g250->Write();
  /*
  g300->Write();
  g400->Write();
  g500->Write();
  g600->Write();
  g700->Write();
  g800->Write();
  g900->Write();
  g1000->Write();
  */
  fout.Close();
  return 1;
}
